#include <kuai/misc/strtool.hpp>
#include <kuai/misc/sqlite3.hpp>
#include <kuai/misc/error.hpp>

#include <kuai/mol/mol.hpp>
#include <kuai/mol/inchi.hpp>
#include <kuai/db/trajdb.hpp>
#include <kuai/mol/molio.hpp>



namespace kuai {

	TrajectoryJob::TrajectoryJob(int mode)
		: DatabaseJob(mode)
	{ }

	void TrajectoryJob::doJob(Sqlite3DatabasePointer pdb, std::ostream& out_log, std::ostream& err_log, StringMap& par) {
		String repository = getParameter(par, "repository", "");
		if (repository.empty()) {
			throw error("Missing Parameter [repository]");
		}
		else {
			pdb->commit();
			String sql = "ATTACH DATABASE \'" + repository + "\' AS repository;";
			pdb->execute(sql.c_str());
			pdb->begin();
		}
	}

	static const char* SQL_SET_CREATE_ANALYSIS_TABLE[] = {
		"CREATE TABLE IF NOT EXISTS repository.moltypes (\n"
			"id INTEGER PRIMARY KEY AUTOINCREMENT,\n"
			"natoms INTEGER,\n"
			"nbonds INTEGER,\n"
			"formula CHAR(50),\n"
			"weight DOUBLE,\n"
			"inchi_key CHAR(50) UNIQUE,\n"
			"data TEXT DEFAULT ''\n"
		");\n",
		"CREATE TABLE molecules (\n"
			"id INTEGER PRIMARY KEY AUTOINCREMENT,\n"
			"moltype INTEGER REFERENCES moltypes (id),\n"
			"atom_hash CHAR(50),\n"
			"bond_hash CHAR(50),\n"
			"data TEXT DEFAULT \'\',\n"
			"UNIQUE (atom_hash, bond_hash)\n"
		");\n",
		"CREATE TABLE mol_in_frame (\n"
			"frame INTEGER REFERENCES frames (id),\n"
			"molecule INTEGER REFERENCES molecules (id),\n"
			"PRIMARY KEY(frame, molecule)\n"
		");\n",
		"CREATE TABLE atom_in_mol (\n"
			"atom INTEGER,\n"
			"molecule INTEGER REFERENCES molecules (id),\n"
			"PRIMARY KEY(atom, molecule)\n"
		");\n",
		"CREATE TABLE bond_in_mol (\n"
			"atom1 INTEGER,\n"
			"atom2 INTEGER,\n"
			"molecule INTEGER REFERENCES molecules (id),\n"
			"PRIMARY KEY(atom1, atom2, molecule)\n"
		");\n"
		/*
		"CREATE VIEW mol_info_view AS\n"
			"SELECT molecules.id as id, natoms, nbonds, formula, weight, inchi_key, moltype, moltypes.data AS inchi_code, \n"
			"atom_hash, bond_hash, min(frame) as birth, max(frame) as death, max(frame)-min(frame) as life\n"
			"FROM molecules INNER JOIN repository.moltypes ON molecules.moltype = moltypes.id\n"
			"INNER JOIN mol_in_frame ON molecules.id = mol_in_frame.molecule GROUP BY (molecules.id);\n",
		"CREATE VIEW frame_statics AS\n"
			"SELECT frame, moltype, natoms, nbonds, formula, weight, inchi_key, inchi_code, count(*) as hits\n"
			"FROM mol_info_view INNER JOIN mol_in_frame ON mol_info_view.id = mol_in_frame.molecule\n"
			"GROUP BY frame, moltype\n"
			"ORDER BY frame, moltype;\n"
		*/
	};

	static const char* SQL_SELECT_ALL_FRAMES = "SELECT frame, data FROM frames;"; 

	static const char* SQL_FIND_MOL_TYPE = "SELECT id FROM moltypes where inchi_key = ?;"; 
	static const char* SQL_INSERT_MOL_TYPE = "INSERT INTO moltypes (natoms, nbonds, formula, weight, inchi_key, data) VALUES (?, ?, ?, ?, ?, ?);"; 

	static const char* SQL_FIND_MOLECULE = "SELECT id FROM molecules WHERE atom_hash = ? AND bond_hash = ?;"; 
	static const char* SQL_INSERT_MOLECULE = "INSERT INTO molecules (moltype, atom_hash, bond_hash) VALUES (?, ?, ?);"; 

	static const char* SQL_INSERT_ATOM_IN_MOL = "INSERT INTO atom_in_mol (atom, molecule) VALUES (?, ?);" ;
	static const char* SQL_INSERT_BOND_IN_MOL = "INSERT INTO bond_in_mol (atom1, atom2, molecule) VALUES (?, ?, ?);";

	static const char* SQL_INSERT_MOL_IN_FRAME = "INSERT INTO mol_in_frame (frame, molecule) VALUES (?, ?);";


	IMPLEMENT_JOB_CREATOR(ParseTrajectoryJob, parse_trajec);

	ParseTrajectoryJob::ParseTrajectoryJob() 
	{
		tor = 0.2f;
	} 

	namespace {

		const char* next_line(const char* p) {
			for (;; ++p) {
				if (*p == NULL) {
					return NULL;
				}
				else if (*p == '\n') {
					return p;
				}
			}
		}

		void parse_mol(const char* value, std::map<Index, RealNumber>& cr, RealNumber tor, PhaseSpace& space) {
			assert (value != NULL && *value != NULL);

			const char* p = next_line(value+1);
			size_t nAtoms = lexical_cast<size_t>(boost::algorithm::trim_copy(String(value, p)));

			std::vector<XYZ> coords;
			std::vector<Index> atoms;
			coords.reserve(nAtoms);
			atoms.reserve(nAtoms);

			value = p;
			p = next_line(value+1);
			PBC* pbc;
			PBC  pbc0;
			try {
				pbc0 = PBC::parse(String(value, p));
				pbc = &pbc0;
			}
			catch (boost::bad_lexical_cast&) {
				pbc = NULL;
			};
			
			while (*value != NULL) {
				value = p;
				p = next_line(value+1);
				if (p == NULL) {
					break;
				}
				else {
					String line(value, p);
					StringArray tokens = split(line);
					if (tokens.size() >= 4) {
						try {
							if (tokens[0].size() == 2) {
								tokens[0][1] = tolower(tokens[0][1]);
							}
							Index e = Element::getElement(tokens[0]).number;
							atoms.push_back(e);
							XYZ v;
							v.x = lexical_cast<RealNumber>(tokens[1]);
							v.y = lexical_cast<RealNumber>(tokens[2]);
							v.z = lexical_cast<RealNumber>(tokens[3]);
							coords.push_back(v);
							continue;
						}
						catch (boost::bad_lexical_cast&) 
						{ }
					}
					throw BadLexicalCast(line, typeid(XYZ));
				}
			}

			IndexArray bonds;
			guessBondLink(atoms.size(), &atoms[0], &coords[0], pbc, bonds, cr, tor);

			std::vector<Index> ct;
			ct.reserve(2 + atoms.size() + bonds.size());
			ct.push_back(atoms.size());
			ct.push_back(bonds.size()/3);
			std::copy(atoms.begin(), atoms.end(), std::back_inserter(ct));
			std::copy(bonds.begin(), bonds.end(), std::back_inserter(ct));

			space[MOL_INFO_ATOM_COORD_3D] = SimObjPtr(new PropertyXYZ(coords, true));
			space[MOL_INFO_CT] = SimObjPtr(new PropertyConnectTable(&ct[0]));
			space[MOL_INFO_PBC] = SimObjPtr(new PropertyPBC(*pbc));
		}
	}

	void ParseTrajectoryJob::createDB() {
		for (size_t i = 0; i < arrayLength(SQL_SET_CREATE_ANALYSIS_TABLE); ++i) {
			pdb->execute(SQL_SET_CREATE_ANALYSIS_TABLE[i]);
		}
		
		int result = pdb->execute(SQL_SELECT_ALL_FRAMES, 
				(Sqlite3Database::ExecCallBack)parseFrame, (void*)this);

		if (result != SQLITE_OK) {
			throw error(errorMessage);
		}
	}


	int ParseTrajectoryJob::parseFrame(ParseTrajectoryJob* pjob, int ncol, char** values, char**) {
		try {
			size_t frame = lexical_cast<size_t>(values[0]);

			Time t1 = now();
			pjob->saveFrame(frame, String(values[1]));
			*(pjob->poutput) << "    Frame " << std::setw(12) << frame << ": " << now() - t1 << std::endl;

			return SQLITE_OK;
		}
		catch (std::exception& error) {
			pjob->errorMessage = error.what();
		}
		catch (...) {
			pjob->errorMessage = SZ_UNKNOWN_ERROR;
		}
		return SQLITE_ABORT;

	}

	void ParseTrajectoryJob::saveFrame(size_t frame, const String& data)
	{
		PhaseSpace space;
		parse_mol(data.c_str(), elementRadii, tor, space);
		ConnectTable& ct = static_cast<PropertyConnectTable&>(*(space.find(MOL_INFO_CT)->second)).x_ct();
		const XYZ* p = static_cast<PropertyXYZ&>(*(space.find(MOL_INFO_ATOM_COORD_3D)->second)).c_ptr();
		const PBC* pbc = NULL;
		if (space.find(MOL_INFO_PBC) != space.end()) {
			pbc = static_cast<PropertyPBC&>(*(space.find(MOL_INFO_PBC)->second)).c_ptr();
		}

		std::vector<boost::shared_ptr<ConnectTable::SubMolType> > subs = getConnectedParts(ct);
		for (std::vector<boost::shared_ptr<ConnectTable::SubMolType> >::const_iterator
			it = subs.begin(); it != subs.end(); ++it)
		{
			Index nAtoms = (*it)->countAtoms(), nBonds = (*it)->countBonds();
			std::vector<Index> atoms; atoms.reserve(nAtoms);
			for (size_t i = 0; i < nAtoms; ++i) {
				atoms.push_back(ct.indexOf((*it)->getAtom(i)));									
			}
			std::sort(atoms.begin(), atoms.end());
			std::vector<std::pair<Index, Index> > bonds; bonds.reserve(nBonds);
			for (size_t i = 0; i < nBonds; ++i) {
				const Bond& bondI = (*it)->getBond(i);
				Index atom1 = ct.indexOf((*it)->getAtom(bondI, 0));
				Index atom2 = ct.indexOf((*it)->getAtom(bondI, 1));
				if (atom1 > atom2) {
					std::swap(atom1, atom2);
				}
				bonds.push_back(std::make_pair(atom1, atom2));
			}
			std::sort(bonds.begin(), bonds.end());
			int molid = saveMolecule(p, **it, ct, atoms, bonds, pbc);

			Sqlite3StatementPointer sqlInsertMolInFrame = pdb->prepare(SQL_INSERT_MOL_IN_FRAME);
			sqlInsertMolInFrame->bind(1, int(frame));
			sqlInsertMolInFrame->bind(2, molid);
			pdb->insert(*(sqlInsertMolInFrame));
		}
	}

	
	void ParseTrajectoryJob::doJob(Sqlite3DatabasePointer pdb, std::ostream& out_log, std::ostream& err_log, StringMap& par) {
		poutput = &out_log;
		perror = &err_log;

		TrajectoryJob::doJob(pdb, out_log, err_log, par);

		molfolder = FileName(getParameter(par, "molfolder", ""));

		tor = getParameter(par, "tor", 0.2f);
		for (StringMap::const_iterator
			it = par.begin(); it != par.end();)
		{
			String symbol = it->first;
			assert (!symbol.empty());
			symbol[0] = ::toupper(symbol[0]);
			const Element& e = Element::getElement(symbol);
			if (e.number != 0) {
				elementRadii[e.number] = lexical_cast<RealNumber>(it->second);
				par.erase(it++);
			}
			else {
				++it;
			}
		}
		createDB();
	}

	int ParseTrajectoryJob::saveMolType(Index nAtoms, Index nBonds, const String& inchi, RealNumber weight, String& formula) {
		String key = inchikey(inchi.c_str());
		key.resize(14);
		String::size_type i1 = inchi.find('/');
		assert (i1 != inchi.npos);
		String::size_type i2 = inchi.find('/', i1+1);
		if (i2 == inchi.npos) {
			i2 = inchi.size();
		}
		formula = inchi.substr(i1+1, i2-i1-1);

		try {
			Sqlite3StatementPointer sqlFindMolType = pdb->prepare(SQL_FIND_MOL_TYPE);
			sqlFindMolType->bind(1, key);
			return int(pdb->getSingleInt(*sqlFindMolType));
		}
		catch (ErrorCanNotFindRow&) {
			Sqlite3StatementPointer sqlInsertMolType = pdb->prepare(SQL_INSERT_MOL_TYPE);
			sqlInsertMolType->bind(1, int(nAtoms));
			sqlInsertMolType->bind(2, int(nBonds));
			sqlInsertMolType->bind(3, formula);
			sqlInsertMolType->bind(4, weight);
			sqlInsertMolType->bind(5, key);
			sqlInsertMolType->bind(6, inchi);
			return int(pdb->insert(*sqlInsertMolType));
		}
	}

	int ParseTrajectoryJob::saveMolecule(const XYZ coords[], SubMolecule<ConnectTable>& mol, const ConnectTable& ct,
		const std::vector<Index>& atoms, const std::vector<std::pair<Index, Index> >& bonds, const PBC* pbc) 
	{
		String md5atoms = md5hex(reinterpret_cast<const Byte*>(&atoms[0]), sizeof(Index)*atoms.size());
		String md5bonds;
		if (!bonds.empty()) {
			md5bonds = md5hex(reinterpret_cast<const Byte*>(&bonds[0]), sizeof(std::pair<Index, Index>)*bonds.size());
		}

		try {
			Sqlite3StatementPointer sqlFindMolecule = pdb->prepare(SQL_FIND_MOLECULE);
			sqlFindMolecule->bind(1, md5atoms);
			sqlFindMolecule->bind(2, md5bonds);
			return int(pdb->getSingleInt(*sqlFindMolecule));
		}
		catch (ErrorCanNotFindRow&) {
			String inchicode[4];
			inchi(mol, "/DoNotAddH /SNon", inchicode);
			RealNumber weight = 0;

			for (size_t i = 0; i < atoms.size(); ++i) {
				weight += Element::getElement(mol.getAtom(i).type).weight;
			}
			String formula;
			int moltype = saveMolType(atoms.size(), bonds.size(), inchicode[0], weight, formula);
			String name = str(moltype) + "." + formula + ".mol";
			if (molfolder.is_dir()) {
				FileName file = molfolder / name;
				saveMolFile(file, coords, mol, ct, pbc);
			}

			int result; 
			{
				Sqlite3StatementPointer sqlInsertMolecule = pdb->prepare(SQL_INSERT_MOLECULE);
				sqlInsertMolecule->bind(1, moltype);
				sqlInsertMolecule->bind(2, md5atoms);
				sqlInsertMolecule->bind(3, md5bonds);
				
				result = int(pdb->insert(*sqlInsertMolecule));
			}

			Sqlite3StatementPointer sqlInsertAtomInMol = pdb->prepare(SQL_INSERT_ATOM_IN_MOL);
			sqlInsertAtomInMol->bind(2, result);
			for (std::vector<Index>::const_iterator
				i = atoms.begin(); i != atoms.end(); ++i)
			{
				sqlInsertAtomInMol->reset();
				sqlInsertAtomInMol->bind(1, int((*i) + 1));
				pdb->insert(*sqlInsertAtomInMol);
			}

			Sqlite3StatementPointer sqlInsertBondInMol = pdb->prepare(SQL_INSERT_BOND_IN_MOL);
			sqlInsertBondInMol->bind(3, result);
			for (std::vector<std::pair<Index, Index> >::const_iterator
				i = bonds.begin(); i != bonds.end(); ++i)
			{
				sqlInsertBondInMol->reset();
				sqlInsertBondInMol->bind(1, int(i->first + 1));
				sqlInsertBondInMol->bind(2, int(i->second + 1));
				pdb->insert(*sqlInsertBondInMol);
			}

			return result;
		}
	
	}

	void ParseTrajectoryJob::saveMolFile(const FileName& file, const XYZ coords[], const SubMolecule<ConnectTable>& mol, const ConnectTable& ct, const PBC* pbc) {
		std::vector<Index> data; data.reserve(2 + mol.countAtoms() + 3 * mol.countBonds());
		std::vector<XYZ> molcoord; molcoord.reserve(mol.countAtoms());

		data.push_back(mol.countAtoms());
		data.push_back(mol.countBonds());
		for (size_t i = 0; i < mol.countAtoms(); ++i) {
			const Atom& atomI = mol.getAtom(i);
			data.push_back(mol.getAtom(i).type);
			Index aI = ct.indexOf(atomI);
			molcoord.push_back(coords[aI]);
		}
		for (size_t i = 0; i < mol.countBonds(); ++i) {
			const Bond& bondI = mol.getBond(i);
			Index a1 = mol.indexOf(mol.getAtom(bondI, 0));
			Index a2 = mol.indexOf(mol.getAtom(bondI, 1));
			data.push_back(bondI.order);
			data.push_back(a1);
			data.push_back(a2);
		}

		PhaseSpace space;
		space[MOL_INFO_CT] = SimObjPtr(new PropertyConnectTable(&data[0]));
		space[MOL_INFO_ATOM_COORD_3D] = SimObjPtr(new PropertyXYZ(molcoord, true));

		if (pbc) {
			pbc->pack_mol(space);
		}

		MolWriterManager::getInstance().write(file, space);	
	};

	JOB_HELP_FUNCTION(ParseTrajectoryJob, output) {
		output << 
			"This function trace the input traject database and Create the ct \n"
			"Parameters:\n"
			"    database=databasename\n"
			"        The name of database\n"
			"    tablename=output_file.txt\n"
			"        tablename=output_file.txt\n"
		;
	}

	static const Char* SQL_CREATE_MOLLIFE_TABLE = 
		"CREATE TABLE mol_life (\n"
			"molecule INTEGER REFERENCES molecules (id),\n"
			"birth INTEGER REFERENCES frames (frame),\n"
			"death INTEGER,\n"
			"PRIMARY KEY(molecule, birth)\n"
		");"
		"CREATE TABLE bond_life (\n"
			"atom1 INTEGER,\n"
			"atom2 INTEGER,\n"
			"birth INTEGER REFERENCES frames (frame),\n"
			"death INTEGER,\n"
			"PRIMARY KEY(atom1, atom2, birth)\n"
		");\n";

	static const Char* SQL_SELECT_TOP2_FRAME = "SELECT frame FROM frames ORDER BY frame LIMIT 2;";

	static const Char* SQL_SELECT_MAX_MOLID = "SELECT max(id) FROM molecules;";
	static const Char* SQL_TRACE_MOL_IN_FRAME = "SELECT frame FROM mol_in_frame WHERE molecule = ? ORDER BY frame;";
	static const Char* SQL_INSERT_MOL_LIFE = "INSERT INTO mol_life (molecule, birth, death) VALUES (?, ?, ?);";
	static const Char* SQL_SELECT_ALL_BOND = "SELECT DISTINCT atom1, atom2 FROM bond_in_mol;";

	static const Char* SQL_TRACE_BOND_IN_FRAME = "SELECT frame FROM mol_in_frame\n"
		"WHERE molecule IN (SELECT molecule FROM bond_in_mol WHERE atom1 = ? AND atom2 = ?)\n"
		"ORDER BY frame";

	static const Char* SQL_INSERT_BOND_LIFE = "INSERT INTO bond_life (atom1, atom2, birth, death) VALUES (?, ?, ?, ?);";

	static int getTimeGap(Sqlite3DatabasePointer pdb) {
		std::vector<StringArray> data;
		pdb->getTable(SQL_SELECT_TOP2_FRAME, data);
		assert (data.size() == 3);
		int i1 = lexical_cast<int>(data[1][0]);
		int i2 = lexical_cast<int>(data[2][0]);
		return i2-i1;
	}

	IMPLEMENT_JOB_CREATOR(TraceMolLifeJob, trace_mollife);

	TraceMolLifeJob::TraceMolLifeJob() 
	{ }
		
	void TraceMolLifeJob::doJob(Sqlite3DatabasePointer pdb, std::ostream& out_log, std::ostream& err_log, StringMap& par) 
	{
		size_t gap = getParameter(par, "gap", size_t(0));
		if (gap == 0) {
			gap = getTimeGap(pdb);
		}
		pdb->execute(SQL_CREATE_MOLLIFE_TABLE);
		saveMolLife(pdb, out_log, gap);
		saveBondLife(pdb, out_log, gap);
	}

	

	void TraceMolLifeJob::saveMolLife(Sqlite3DatabasePointer pdb, size_t i, size_t gap) 
	{
		std::vector<std::pair<size_t, size_t> > life;
		life.reserve(1024);
		Sqlite3StatementPointer pq = pdb->prepare(SQL_TRACE_MOL_IN_FRAME);
		pq->bind(1, sqlite3_int64(i));
		int tag = pq->step();
		if (tag == SQLITE_ROW) {
			size_t finish = size_t(pq->getInt64(0));
			size_t start;
			for (;;) {
				start = finish;
				size_t v2 = 0;
				for (;;) {
					tag = pq->step();
					if (tag == SQLITE_ROW) {
						v2 = size_t(pq->getInt64(0));
						if (v2 != finish + gap) {
							break;
						}
						else {
							finish = v2;
						}						
					}
					else {
						break;
					}
				}
				assert (start <= finish);
				life.push_back(std::make_pair(start, finish));
				if (tag == SQLITE_DONE) {
					break;
				}
				else {
					start = finish = v2;
				}
			}
		}
		pq = pdb->prepare(SQL_INSERT_MOL_LIFE);
		pq->bind(1, sqlite3_int64(i));
		for (size_t j = 0; j < life.size(); ++j) {
			pq->bind(2, sqlite3_int64(life[j].first));
			pq->bind(3, sqlite3_int64(life[j].second));
			pq->step();
			pq->reset();
		}
	}

	void TraceMolLifeJob::saveMolLife(Sqlite3DatabasePointer pdb, std::ostream& out_log, size_t gap) 
	{
		size_t nMols = size_t(pdb->getSingleInt(SQL_SELECT_MAX_MOLID));
		for (size_t i = 1; i <= nMols; ++i) {
			Time t1 = now();
			saveMolLife(pdb, i, gap);
			out_log << "    MolLife " << std::setw(12) << i << ": " << now() - t1 << std::endl;
		}
	
	}
	void TraceMolLifeJob::saveBondLife(Sqlite3DatabasePointer pdb, std::ostream& out_log, size_t gap) 
	{
		std::vector<std::pair<Index, Index> > bonds; bonds.reserve(2048);
		{
			Sqlite3StatementPointer pq = pdb->prepare(SQL_SELECT_ALL_BOND);
			while (pq->step() != SQLITE_DONE) {
				bonds.push_back(std::make_pair(Index(pq->getInt64(0)), Index(pq->getInt64(1))));
			}
		}
		for (size_t i = 0; i < bonds.size(); ++i) {
			Time t1 = now();
			saveBondLife(pdb, bonds[i].first, bonds[i].second, gap);
			out_log << "    BondLife " << std::setw(6) << bonds[i].first 
					<< std::setw(6) << bonds[i].second << ": " << now() - t1 << std::endl;
		}
	}

	void TraceMolLifeJob::saveBondLife(Sqlite3DatabasePointer pdb, Index a1, Index a2, size_t gap) {
		std::vector<std::pair<size_t, size_t> > life;
		life.reserve(1024*2);
		Sqlite3StatementPointer pq = pdb->prepare(SQL_TRACE_BOND_IN_FRAME);
		pq->bind(1, sqlite3_int64(a1));
		pq->bind(2, sqlite3_int64(a2));
		int tag = pq->step();
		if (tag == SQLITE_ROW) {
			size_t finish = size_t(pq->getInt64(0));
			size_t start;
			for (;;) {
				start = finish;
				size_t v2 = 0;
				for (;;) {
					tag = pq->step();
					if (tag == SQLITE_ROW) {
						v2 = size_t(pq->getInt64(0));
						if (v2 != finish + gap) {
							break;
						}
						else {
							finish = v2;
						}						
					}
					else {
						break;
					}
				}
				if (start != finish) {
					assert (start < finish);
					life.push_back(std::make_pair(start, finish));
				}
				if (tag == SQLITE_DONE) {
					break;
				}
				else {
					start = finish = v2;
				}
			}
		}
		pq = pdb->prepare(SQL_INSERT_BOND_LIFE);
		pq->bind(1, sqlite3_int64(a1));
		pq->bind(2, sqlite3_int64(a2));
		for (size_t j = 0; j < life.size(); ++j) {
			pq->bind(3, sqlite3_int64(life[j].first));
			pq->bind(4, sqlite3_int64(life[j].second));
			pq->step();
			pq->reset();
		}
	}

	JOB_HELP_FUNCTION(TraceMolLifeJob, output) {
		output << 
			"This function trace the input traject database and Create the ct \n"
			"Parameters:\n";
		DatabaseJob::help(output);
		output << 
			"    database=databasename\n"
			"        The name of database\n"
			"    tablename=output_file.txt\n"
			"        tablename=output_file.txt\n"
		;
	}



	IMPLEMENT_JOB_CREATOR(TraceReactionJob, trace_reaction);

	static const Char* SQL_CREATE_RXN_TABLE[] = {
		"CREATE TABLE IF NOT EXISTS repository.rxntypes (\n"
			"id INTEGER PRIMARY KEY AUTOINCREMENT,\n"
			"reactants INTEGER,\n"
			"products  INTEGER,\n"
			"formula TEXT,\n"
			"inchi_key CHAR(50) UNIQUE,\n"
			"data TEXT\n"
		");",
		"CREATE TABLE IF NOT EXISTS repository.moltype_in_rxntype (\n"
			"rxntype INTEGER REFERENCES rxntypes(id), \n"
			"moltype INTEGER REFERENCES moltypes(id), \n"
			"character CHAR(1),\n"
			"PRIMARY KEY (rxntype, moltype, character)\n"
		");",
		"CREATE TABLE reactions (\n"
			"id INTEGER PRIMARY KEY AUTOINCREMENT,\n"
			"rxntype INTEGER REFERENCES rxntypes (id),\n"
			"frame INTEGER REFERENCES frames (frame)\n"
		");",
		"CREATE TABLE mol_in_rxn (\n"
			"reaction INTEGER REFERENCES reactions(id), \n"
			"molecule INTEGER REFERENCES molecules(id), \n"
			"character CHAR(1),\n"
			"PRIMARY KEY (reaction, molecule, character)\n"
		");",
	};

	static const Char* SQL_GET_FIRST_FRAME = "SELECT frame FROM frames ORDER BY frame LIMITS 1;";

	TraceReactionJob::TraceReactionJob()
	{ }
	
	void TraceReactionJob::doJob(Sqlite3DatabasePointer pdb, std::ostream& out_log, std::ostream& err_log, StringMap& par)
	{
		TrajectoryJob::doJob(pdb, out_log, err_log, par);
		plog = &out_log;
		size_t gap = getParameter(par, "gap", size_t(0));
		if (gap == 0) {
			gap = getTimeGap(pdb);
		}
		loadMolecule();

		for (size_t i = 0; i < arrayLength(SQL_CREATE_RXN_TABLE); ++i) {
			pdb->execute(SQL_CREATE_RXN_TABLE[i]);
		}
		std::vector<int> rxn_frames;
		pdb->getIntArray("SELECT DISTINCT birth FROM mol_life ORDER BY birth", 0, rxn_frames);
		for (size_t i = 1; i < rxn_frames.size(); ++i) {
			int birth = rxn_frames[i];
			int death = birth - gap;
			saveRxn(birth, death);
		}

	}

	void TraceReactionJob::saveRxn(int birth, int death) {
		std::vector<int> reactants, products;
		pdb->getIntArray("SELECT molecule FROM mol_life WHERE death = " + str(death), 0, reactants);
		pdb->getIntArray("SELECT molecule FROM mol_life WHERE birth = " + str(birth), 0, products);

		if (reactants.size() == 1 || products.size() == 1) {
			saveRxn(reactants, products, birth);
		}
		else {
			splitRxn(reactants, products, birth);
		}
	}

	void TraceReactionJob::saveRxn(const std::vector<int>& reactants, const std::vector<int>& products, int frame) 
	{
		StringArray kR, kP;
		for (size_t i = 0; i < reactants.size(); ++i) {
			kR.push_back(mols[reactants[i]].key);
		}
		for (size_t i = 0; i < products.size(); ++i) {
			kP.push_back(mols[products[i]].key);
		}
		std::sort(kR.begin(), kR.end());
		std::sort(kP.begin(), kP.end());
		String coder;
		for (size_t i = 0; i < kR.size(); ++i) {
			coder += kR[i];
		}
		coder += "=";
		for (size_t i = 0; i < kP.size(); ++i) {
			coder += kP[i];
		}

		String key = md5hex(coder);

		int rxntype;
		try {
			rxntype = int(pdb->getSingleInt("SELECT id FROM rxntypes WHERE inchi_key = \'" + key + "\'"));
		}
		catch (ErrorCanNotFindRow&) {
			rxntype = saveRxnType(reactants, products, key);

			coder = "";
			for (size_t i = 0; i < kP.size(); ++i) {
				coder += kP[i];
			}
			coder += "=";
			for (size_t i = 0; i < kR.size(); ++i) {
				coder += kR[i];
			}
			String key2 = md5hex(coder);
			if (key2 == key) {
				key2 = "<" + key2 + ">";
			}
			saveRxnType(products, reactants, key2);
		}

		int rxnid = int(pdb->insert("INSERT INTO reactions (rxntype, frame) VALUES (" + str(rxntype) + ", " + str(frame) + ")"));

		Sqlite3StatementPointer p = pdb->prepare("INSERT INTO mol_in_rxn (reaction, molecule, character) VALUES (?, ?, ?)");
		p->bind(1, rxnid);
		p->bind(3, "R");
		for (size_t i = 0; i < reactants.size(); ++i) {
			p->bind(2, reactants[i]);
			pdb->insert(*p);
			p->reset();
		}
		p->bind(3, "P");
		for (size_t i = 0; i < products.size(); ++i) {
			p->bind(2, products[i]);
			pdb->insert(*p);
			p->reset();
		}
		*plog << "Save Reaction " << rxnid << ": " << key << " of " << rxntype << " @ " << frame << std::endl;
	}

	int TraceReactionJob::saveRxnType(const std::vector<int>& reactants, const std::vector<int>& products, const String& key)
	{
		std::map<int, int> rc, pc;
		std::map<int, String> formula;
		for (size_t i = 0; i < reactants.size(); ++i) {
			int itype = mols[reactants[i]].type;
			formula[itype] = mols[reactants[i]].formula;
			++rc[itype];
		}
		for (size_t i = 0; i < products.size(); ++i) {
			int itype = mols[products[i]].type;
			formula[itype] = mols[products[i]].formula;
			++pc[itype];
		}
		String f1, f2;
		for (std::map<int, int>::const_iterator
			it = rc.begin(); it != rc.end(); ++it)
		{
			if (it->second > 1) {
				f1 += str(it->second);
				f2 += str(it->second);
			}
			f1 += "[" + formula[it->first] + "]+";
			f2 += "[" + str(it->first) + "]+";
		}
		assert (!f1.empty() && !f2.empty());
		f1.resize(f1.size()-1);
		f2.resize(f2.size()-1);
		f1 += " = ";
		f2 += " = ";

		for (std::map<int, int>::const_iterator
			it = pc.begin(); it != pc.end(); ++it)
		{
			if (it->second > 1) {
				f1 += str(it->second);
				f2 += str(it->second);
			}
			f1 += "[" + formula[it->first] + "]+";
			f2 += "[" + str(it->first) + "]+";
		}
		f1.resize(f1.size()-1);
		f2.resize(f2.size()-1);

		Sqlite3StatementPointer p1 = pdb->prepare("INSERT INTO rxntypes (reactants, products, formula, inchi_key, data) VALUES(?, ?, ?, ?, ?);");
		p1->bind(1, int(reactants.size()));
		p1->bind(2, int(products.size()));
		p1->bind(3, f1);
		p1->bind(4, key);
		p1->bind(5, f2);
		int result = int(pdb->insert(*p1));
		*plog << "Save # Type # " << result << ": " << f1 << "  | " << f2 << std::endl;

		p1 = pdb->prepare("INSERT INTO moltype_in_rxntype (rxntype, moltype, character) VALUES(?, ?, ?);");
		p1->bind(1, result);
		p1->bind(3, "R");
		for (std::map<int, int>::const_iterator 
			it = rc.begin(); it != rc.end(); ++it)
		{
			p1->bind(2, it->first);
			p1->step();
			p1->reset();
		}


		p1->bind(3, "P");
		for (std::map<int, int>::const_iterator 
			it = pc.begin(); it != pc.end(); ++it)
		{
			p1->bind(2, it->first);
			p1->step();
			p1->reset();
		}
		return result;
	}

	void TraceReactionJob::loadMolecule() {
		std::vector<StringArray> data;
		data.reserve(2048);
		pdb->getTable("SELECT molecules.id as id, inchi_key, formula, moltype FROM molecules INNER JOIN repository.moltypes ON molecules.moltype = moltypes.id ORDER BY id", data);

		mols.resize(data.size());
		for (size_t i = 1; i < data.size(); ++i) {
			assert (str(i) == data[i][0]);
			mols[i].key = data[i][1];
			mols[i].formula = data[i][2];
			mols[i].type = lexical_cast<int>(data[i][3]);
		}
		Sqlite3StatementPointer p = pdb->prepare("SELECT molecule, atom FROM atom_in_mol");
		while (p->step() == SQLITE_ROW) {
			int mol = int(p->getInt64(0));
			int atom = int(p->getInt64(1));
			assert (1 <= mol && mol <= int(mols.size()));
			mols[mol].atoms.push_back(atom);
		}
	}

	void TraceReactionJob::splitRxn(const std::vector<int>& reactants, const std::vector<int>& products, int frame) 
	{
		std::vector<bool> searched(reactants.size(), false);
		for (size_t i = 0; i < searched.size(); ++i) {
			if (searched[i] == false) {
				std::set<int> r, p;
				std::set<int> ratom, patom;
				searched[i] = true;
				int rmol1 = reactants[i];
				r.insert(rmol1);
				ratom.insert(mols[rmol1].atoms.begin(), mols[rmol1].atoms.end());
				while (ratom != patom) {
					for (size_t j = 0; j < products.size(); ++j) {
						int pmolJ = products[j];
						if (hasSameAtom(ratom, mols[pmolJ].atoms)) {
							patom.insert(mols[pmolJ].atoms.begin(), mols[pmolJ].atoms.end());
							p.insert(pmolJ);
						}
					}
					for (size_t j = i+1; j < reactants.size(); ++j) {
						if (!searched[j]) {
							int rmolJ = reactants[j];
							if (hasSameAtom(patom, mols[rmolJ].atoms)) {
								searched[j] = true;
								ratom.insert(mols[rmolJ].atoms.begin(), mols[rmolJ].atoms.end());
								r.insert(rmolJ);
							}
						}
					}
				}
				saveRxn(std::vector<int>(r.begin(), r.end()), std::vector<int>(p.begin(), p.end()), frame);
			}
		}
	}

	bool TraceReactionJob::hasSameAtom(const std::set<int>& atomI, const std::vector<int>& atomJ) {
		for (size_t j = 0; j < atomJ.size(); ++j) {
			if (atomI.find(atomJ[j]) != atomI.end()) {
				return true;
			}
		}
		return false;
	}


	JOB_HELP_FUNCTION(TraceReactionJob, output) {
		output << 
			"This function trace the input traject database and Create the ct \n"
			"Parameters:\n";
		DatabaseJob::help(output);
		output << ""
		;	
	}

	IMPLEMENT_JOB_CREATOR(DumpTrajectoryMoleculeJob, dump_traj_mol);

	DumpTrajectoryMoleculeJob::DumpTrajectoryMoleculeJob()
		: TrajectoryJob(SQLITE_OPEN_READONLY)
	{ }

	void DumpTrajectoryMoleculeJob::doJob(Sqlite3DatabasePointer pdb, std::ostream& out_log, std::ostream& err_log, StringMap& par) {
		// TODO:

	};

	JOB_HELP_FUNCTION(DumpTrajectoryMoleculeJob, output) {
		output << 
			"This function dump molecule in trajectory database\n"
			"Parameters:\n";
		DatabaseJob::help(output);
		output << ""
		;	
	}

	IMPLEMENT_JOB_CREATOR(ParseReactionGraph, reaction_graph);

	void ParseReactionGraph::doJob(std::ostream& osOutput, std::ostream& osError, StringMap& par) {
		String file = getParameter(par, "reaction");
		std::ifstream stream(file.c_str());

		String line;
		std::map<Pair, int> weights;
		std::vector<Pair> edges;
		int nMol = 0;
		while (getline(stream, line)) {
			StringArray tokens = split(line, "|");
			int w = lexical_cast<int>(tokens[3]);
			std::pair<String, String> rxn = splitPair(tokens[2], '=');
			StringArray r = split(rxn.first, "+");
			StringArray p = split(rxn.second, "+");

			std::vector<int> factorR;
			std::vector<int> molR;
			for (size_t i = 0; i < r.size(); ++i) {
				r[i] = strip(r[i]);
				int factor = 1;
				if (r[i][0] != '[') {
					size_t ii = r[i].find('[');
					factor = lexical_cast<size_t>(r[i].substr(0, ii));
					r[i] = r[i].substr(ii);
				}
				r[i] = r[i].substr(1, r[i].size()-2);
				int ri = lexical_cast<int>(r[i]);
				factorR.push_back(factor);
				molR.push_back(ri);
				if (ri > nMol) {
					nMol = ri;
				}
			}

			std::vector<int> factorP;
			std::vector<int> molP;
			for (size_t j = 0; j < p.size(); ++j) {
				p[j] = strip(p[j]);
				int factor = 1;
				if (p[j][0] != '[') {
					size_t jj = p[j].find('[');
					factor = lexical_cast<size_t>(p[j].substr(0, jj));
					p[j] = p[j].substr(jj);
				}
				p[j] = p[j].substr(1, p[j].size()-2);
				int pj = lexical_cast<int>(p[j]);
				factorP.push_back(factor);
				molP.push_back(pj);

				if (pj > nMol) {
					nMol = pj;
				}
			}

			for (size_t i = 0; i < molR.size(); ++i) {
				for (size_t j = 0; j < molP.size(); ++j) {
					if (molR[i] != molP[j]) {
						Pair e = Pair(molR[i], molP[j]);
						edges.push_back(e);
						weights[e] = factorR[i] * factorP[j] * w;
					}
				}
			}
		}
				/*
		Graph g(edges.begin(), edges.end(), nMol+1);
		boost::graph_traits<Graph>::vertex_iterator vi, vi_end, next;
		boost::tie(vi, vi_end) = vertices(g);
		bool changed=true;
		while (changed) {
			changed=false;
			for (next = vi; vi != vi_end; vi = next) {
				++next;
				if (boost::in_degree(*vi, g) == 0) {
					remove_vertex(*vi, g);
					changed=true;
				}
			}
		}


		boost::graph_traits<Graph>::edge_iterator ei, ei_end;
		boost::tie(ei, ei_end) = boost::edges(g);
		for (; ei != ei_end; ++ei) {
			osOutput << source(*ei, g) << ", " << target(*ei, g) << std::endl;
		}
		*/




		// parseData(g, osOutput);
	}

	void ParseReactionGraph::parseData(ParseReactionGraph::Graph& g, std::ostream& osOutput) {
		typedef std::vector< Vertex > container;
		container c;
		topological_sort(g, std::back_inserter(c));

		osOutput << "A topological ordering: ";
		for ( container::reverse_iterator ii=c.rbegin(); ii!=c.rend(); ++ii) {
			osOutput << *ii << " ";
		}
		osOutput << std::endl;	
	}

	JOB_HELP_FUNCTION(ParseReactionGraph, output) {
		output << 
			"This function try to create the graph of reactions\n"
			"Parameters:\n";
	}

}


